function B=MagneticField_ESL(x_d,y_d,z_d,l,Jedw,ae)

r_d=sqrt(x_d.^2+y_d.^2);

if r_d>0
    cos_theta_d=x_d/r_d;
    sin_theta_d=y_d/r_d;
    
    B_R=dblquad(@(z,theta)MagneticField3D_R_Func...
        (ae,z,theta,r_d,z_d),-l/2,l/2,0,2*pi);
else
    B_R=0;
    cos_theta_d=0;
    sin_theta_d=0;
end

B_Z=dblquad(@(z,theta)MagneticField3D_Z_Func...
    (ae,z,theta,r_d,z_d),-l/2,l/2,0,2*pi);

B=-Jedw*[B_R*cos_theta_d B_R*sin_theta_d B_Z];